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We show that the planetesimal formation due to the gravitational fragmen- 
tation of a dust layer in a protoplanetary disk is possible. The dust density 
distribution in the dust layer would approach the constant Richardson number 
To . distribution due to the dust stirring by the shear instability and dust settling. We 

perform the analysis of the shear instability of dust layer in a protoplanetary disk 



with the constant Richardson number density distribution. Our study revealed 
that this distribution is stable against the shear instability even if the dust den- 
sity at the midplane reaches the critical density of the gravitational instability, 
and the planetesimal formation through the gravitational fragmentation of the 
dust layer can occur even for the dust to gas surface density ratio with the solar 
composition. 

Subject headings: planetary systems: protoplanetary disks — solar system: formation- 
hydrodynamics — instabilities 



1. Introduction 

It is in controversy whether planetesimals are formed in a protoplanetary disk by the 
gravitational instability (Safronov 1969; Goldreich & Ward 1973; Coradini, Feferico, & Magni 
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1981; Sekiya 1983) or by collisional sticking of dust aggregates due to some non-gravitational 
forces (Weidenschilling 1984; Cuzzi, Dobrovolskis, & Champney 1993; Wurm, Blum, & Col- 
well 2001). 

The planetesimal formation through dust sticking due to non-gravitational forces has 
some difficulties. A meter-sized body falls to the inner edge of the disk due to the gas drag 
in much shorter time scale (~ 10'^ yr) (Adachi, Hayashi, & Nakazawa 1976; Weidenschilling 
1977) than the possible life time of a disk (~ 10^ yr). Cuzzi, Dobrovolskis, & Champney 
(1993) have suggested that a body sweeps smaller dust aggregates as the body falls towards 
the central star and can grow to form a km-sized planetesimal. Moreover, Wurm, Blum, & 
Colwell (2001) have observed the growth of a body due to the following mechanism: When 
a dust aggregate impacts a target body, the aggregate fragments into monomers. After the 
monomers bound on the body several times, they finally stick to the body. However, this 
experiment was done with a condition in which the target body was smaller than the gas 
mean free path. Sekiya & Takeda (2003) have indicated that the dust aggregate breaks into 
monomers by the first collision and the monomers do not collide the target body again if 
the body is larger than the gas mean free path, because the monomers are carried away by 
the gas fiow. Thus, the body can no longer grow by sweeping monomers if the body grows 
up to a radius on the order of the mean free path of the gas. For example, the body can not 
grow larger than about 1cm at lAU. 

On the other hand, the formation of planetesimals due to the gravitational instabil- 
ity is attractive since the particular sticking mechanism of dust aggregates is not needed. 
Moreover, by growing from millimeter-sized dust aggregates to a km-sized body within a 
short period, the process avoids the problem that meter-sized bodies fall towards the central 
star by the gas drag. However, it has been considered that the gravitational instability of 
the dust layer is difficult to occur if a disk is turbulent; the dust aggregates are stirred up 
from the midplane, so that the dust layer can not reach the critical density above which the 
gravitational instability of the dust layer arises. 

The local turbulence would occur even though the disk is laminar, that is, the global tur- 
bulence such as the thermal convection does not exist. Weidenschilling (1980) has suggested 
that the turbulence is induced due to the vertical shear in the dust layer. The turbulence 
diffuses the dust aggregates and prevents them from settling towards the midplane. The 
cause of the vertical shear is explained as follows. The gas tends to revolve with a sub- 
Kepler velocity, balancing the gravitation of the central star, the centrifugal force, and the 
gas pressure gradient. On the other hand, dust tends to revolve with the Kepler velocity, 
balancing the two former forces. Thus, supposing the coupling between dust and gas is good. 
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their velocity is determined by 

V4>=[^-{'nPg/p)\vK, (1) 

where 

P = Pg + Pd, (2) 

V = -r{dP/dr)/{2vj,Pg), (3) 

where pg is the spatial gas density, pd is the spatial dust density (i.e. the total mass of dust 
in unit volume of the protoplanetary disk), vk is the Kepler velocity, r is the distance from 
the rotation axis of the disk, and P is the gas pressure. Thus shear occurs in the dust layer, 
depending on the dust to gas density ratio (Weidenschilling 1980; Cuzzi, Dobrovolskis, & 
Champney 1993; Sekiya 1998). 

Cuzzi, Dobrovolskis, & Champney (1993) and Sekiya (1998) have argued that the dust 
diffusion by shear-induced turbulence and the dust settling by the vertical component of the 
gravity of a central star lead to an equilibrium dust distribution. The strength of the shear 
is characterized by the Richardson number J, which is defined by 

J = -g,{l/p){dp/dz){dv^/dz)-\ (4) 

where z is the distance from the midplane, and Qz is the gravitational acceleration in the z- 
direction. In a shear flow with density stratification and the gravitational force, but without 
the effects of the Coriolis and the tidal forces, the Richardson number J is used in the 
criterion for stability. The physical meaning of J is a quarter of the ratio of the potential 
energy consumed by the perturbation to the extra energy gained by the perturbation from 
the unperturbed shear flow. If J > Jc = 1/4, more energy is consumed than is gained by 
the perturbation; such a motion is inhibited by the energy conservation law and thus the 
flow is stable. On the other hand, the flow can be unstable below the critical Richardson 
number J^. Sekiya (1998) derived the equilibrium state as follows: As dust aggregates settle 
towards the midplane, the local Richardson number decreases below the critical value Jc in 
a certain region of the dust layer. Then the turbulence is induced by the shear instability, 
which would stir dust aggregates. As the dust aggregates diffuse, the shear rate \dv^/dz\ 
decreases, the local Richardson number increases, and the shear-induced turbulence ceases 
after J exceeds the critical value Jc- Then the dust settling starts again. Repeat of the 
processes is expected to bring to the constant Richardson number dust density distribution 
(CRNDDD in the followings) with J = Jc for the entire region of the dust layer if dust 
aggregates are sufficiently small and are stirred up by very weak turbulence. 

The dust density can not reach the critical density of the gravitational instability for 
CRNDDD with J = 0.25 if the dust to gas surface density ratio is that calculated from the 
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solar abundance. In contrast, the gravitational instability could occur when the dust to gas 
surface density ratio increases, that is, the dust concentrates or the gas depletes by some 
mechanisms (Sekiya 1998; Youdin & Shu 2002; Youdin & Chiang 2004). 

In order to elucidate the shear instability in more detail, we have performed a series of 
researches. Sekiya & Ishitsu (2000) have performed the linear analysis of the shear instability 
for CRNDDD. Their results have shown that the growth rate of the shear instability is less 
than the Kepler angular frequency. Next we investigated the instability using a different 
unperturbed state (Sekiya & Ishitsu 2001). As dust particles stick each other and settle 
towards the midplane, the density distribution in the dust layer tends to be constant ( 
Watanabe & Yamada 2000). Thus, we employed the hybrid dust density distribution (called 
HDDD in the following) in which the dust layer has a constant density and transition regions 
have sinusoidal density distributions. For this distribution, the growth rate of the shear 
instability is much larger than the Kepler angular frequency when the dust density at the 
midplane is larger than the gas density; the latter is much less than the critical density of 
the gravitational instability. 

However these works did not include the Coriolis force and the tidal force. Ishitsu & 
Sekiya (2002) have investigated the stability of HDDD including the effect of the Coriolis 
force but neglecting the tidal force. Their results have shown that the Coriolis force has 
little effect on the growth rate of the shear instability for HDDD. However, the mode of the 
instability changes from the co-rotation resonance to resonances like the Lindblad resonance. 

Ishitsu & Sekiya (2003) have investigated the shear instability of HDDD including both 
the tidal force and the Coriolis force. That is, they have employed the dust distribution 
expected before the onset of the shear instability. The paper indicated that the tidal force 
stabilizes the shear instability if the growth rate of the shear instability in the case without 
the tidal force is comparable to or less than the Kepler angular velocity. However, the shear 
instability occurs before the dust density reaches the critical density of the gravitational 
instability. Thus, turbulence due to the shear instability would change the dust density 
distribution. As a result, the dust distribution is expected to approach CRNDDD asymptot- 
ically. In this work, we investigate the linear stability of CRNDDD, including the Coriolis 
and tidal forces and the self-gravity of the dust layer. 

In §2, the constant Richardson number distribution is introduced. In §3, the basic 
equations for the linear analysis are derived. In §4, results of the case with the Coriolis force 
only are shown. In §5, results of the case with the Coriolis and tidal forces are shown. The 
results of the case with the self-gravity of the dust layer are also shown. In §6, possibility of 
the planetesimal formation due to the gravitational instability is discussed. 
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2. The Constant Richardson Number Dust Density Distribution 

Here, the constant Richardson number dust density distribution (CRNDDD) derived 
by Sekiya (1998) is reviewed. As written in the previous section, Sekiya (1998) considered 
that the dust density distribution asymptotically approaches CRNDDD with J = J^, and he 
analytically derived it (note that his solution is also applicable to the case without J = J^, 
thus we here assume only that J =constant): 
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where 
and 



u = pgl\pd{z)^ pc^ + g, 



(5) 

(6) 

(7) 



q = AnGpg/Ql. 

Here, G is the gravitational constant, Qk is the Kepler angular velocity, and Um is the value 
of u at the midplane z = 0. The assumption that the gas density pg is constant within the 
dust layer is appropriate if the thickness of the dust layer is much thinner than the gas scale 
height Hg. Then gas density is given by 

Pg = J^g/iV^Hg), (8) 

where the disk is assumed isothermal in the z-direction and S^ is the gas surface density: 



1.7 X 10^9 (r/AU) -3/2 



g cm 



-21 



(9) 



where fg is the ratio of our gas surface density to that of the Hayashi model (Hayashi 1981; 

Hayashi, Nakazawa, & Nakagawa 1985). The dust surface density is given by 

rzd 
Sd = 2 / pddz 



= 2y/jr]rpg |(1 + g) In [(l + g + ^(1 + Q? - ul) /m^J - ^(1 + q)^ - u. 
= 7.lUr/AVr'/\cm-% 
where the half of thickness of the dust layer z^ is given from equation (5) by 

"m-^ln (l + g + V(l + g)2 - u^a) /u, 



Zd/{^r]r) = V(l + g)^ 



(10) 



;ii) 



and fd is the ratio of the dust surface density to the Hayashi model for the region where the 
water vaporizes. In the disk with the gas to dust mass ratio of the solar abundance. 



fd/fg 



4.2 for the region where water condenses 
1 for the region where water vaporizes. 



(12) 
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Sekiya (1998) considered that CRNDDD with J = J^ is reahzed as written in §1. 
He showed that the midplane density for CRNDDD with J = J^ is much lower than the 
critical density of the gravitational instability. He also showed that the the midplane density 
increases if fd/ fg increases by dust concentration and/or gas depletion by some processes. 
Note here that the critical density of the gravitational instability used by Sekiya (1998) 
Pd{0)/Pg = 260/~^(r/AU)~^/'^ is that for the constant dust density distribution (Sekiya 1983), 
although Sekiya (1998) investigated the stability of CRNDDD. Recently, Yamoto & Sekiya 
(2004) calculated the critical density of the gravitational instability for CRNDDD: Pd{0)/Pg = 
750/_"^(r/AU)~^/'^. The explanation for this difference is that the effective thickness of 
the dust layer for CRNDDD is thinner than the thickness for the constant dust density 
distribution if the dust and gas surface density are same. Thus the CRNDDD requires 
the larger critical density for the gravitational instability. Using this new value, the dust 
density at the midplane reaches the critical density if the dust surface density is enhanced 
by fd = 17.15, or the gas surface density is reduced by fg = 0.0286 compared to Hayashi's 
dust surface density at lAU (Hayashi 1981; Hayashi, Nakazawa, & Nakagawa 1985). These 
values are not so different from those obtained by Sekiya (1998), fd = 16.8 and fg = 0.029, 
because the midplane density increases rapidly as fd increases or fg decreases (Yamoto & 
Sekiya 2004). 



3. Basic Equations 

This section gives the basic equations under some assumptions described below; they are 
similar to those written in Ishitsu & Sekiya (2003), but include the self-gravity. We neglect 
the curvature of the cylindrical coordinates (r, 0, z) and use the local Cartesian coordinate 
system which rotates with the Kepler velocity vk{R) aXr = R where R =constant (Goldreich 
& Tremaine 1978). Our coordinates x,y and z denote the radial, the azimuthal and the 
vertical directions of the disk, respectively. That is, x = r — R, y = R[(f) — QK{R)t], and 
z = z, where Qk{R) is the Kepler angular velocity at r = i? and we neglect higher order terms 
of X, y and z. In the following, we denote vk{R) and Qk{R) by vk and Qk, respectively, 
for simplicity. In the Keplerian disk, the flow experiences the tidal force 3^2^-0;. The gas 
can be assumed incompressible because the dust layer treated here is much thinner than 
vertical scale height of the gas disk. The mixture of gas and dust is considered as one fluid 
because we treat dust aggregates whose frictional times are much shorter than the Kepler 
period and the growth time of the shear instability. This assumption is fairly good if dust 
aggregates are smaller than e.g. 1cm at 1 AU. Thus, we obtain the continuity equation, the 
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mass conservation equation, the momentum equations, and the Poisson's equation 

du dv dw 
dx dy dz 

dp dp dp dp ^ ,^ ,, 

dt dx dy dz 
du du du du IdP ^^^^2 ^^r^ <9\1/ .^^. 

— + M— + 1;— + w— = — ^— + scrn^x + 2cnKV - ^— , 15 

dt dx dy dz p dx dx 

dv dv dv dv 1 dP „^^ (9\E' ,^„, 

^ + M— + t;— + w— = — -2CnKU---, (16) 

dt dx dy dz p oy dy 

dw dw dw dw IdP 2 "^^ /-,7\ 

dt dx dy dz p dz ^ dz^ 

d^ ^ ^-4^^- fl8) 

dx"^ dy"^ dz"^ 

where m, v and w are the radial, the azimuthal and the vertical velocities, respectively, and 
\E' is the gravitational potential. If C = 1 and T = 0, the Coriolis force is taken into account, 
but the tidal force is not taken into account. On the other hand, if C = T = 0, neither the 
Coriolis force nor the the tidal force are taken into account. Of course, T = C = 1 for a real 
flow, but we used these parameters to investigate the effects of the tidal and the Coriolis 
forces by comparing the results with T = C = 0, T = and C = 1, and T = C = 1. For 
u = P = "^ = 0, equation (15) reads 

3 
v = — T^kx. (19) 

This expresses the circular Kepler motion in the local coordinate system. In order to elim- 
inate the Keplerian part of the velocity, we introduce the velocity relative to the Keplerian 

motion, 

3. 



v = v^ -T^kx. (20) 



From equations (13) to (18), we have 



du dv dw , , 

dx dy dz 

du du , 3 _^ ,du du 1 dP ^^ d"^ , , 

— +u— + {v- -T^Kx)— + w— = —— + 2CnKV - — , 23 

dt dx 2 dy dz p dx dx 



dv 
dt 


dv 
dx 


+ {v 




dw 


dw 



3_^ ,dv dv IdP , ^ 3_,^ d^ , ,, 

_™,.)- + u>^ = -- g^ - (2C - -T)n..„, - g^. (24) 

.^ , 3^^ .dw dw IdP ^2 .9^ 

a^ 9^ ^ = 4^^- f26) 

dx"^ dy"^ dz^ 

In order to investigate the linear stability, we assume that an unperturbed state is steady 
and uniform in x and y directions: 

l = ± = ± = 0, (27) 

dt dx dy 

except for dPo/dx =constant7^ 0. We also assume that the unperturbed velocity is parallel 
to the y-axis (i.e. the azimuthal direction) 

uo = wo = 0, (28) 

and pg =constant throughout the region of consideration. The unperturbed azimuthal ve- 
locity is determined from equation (1), which is rewritten 

^0 = -pgiivx/po, (29) 

where r] is assumed to be constant in our local Cartesian coordinates. 

From equations (23), (27) and (28), we have 

1 dP 

_ = 2CQkVo. (30) 

Po ox 

Integrating equation (26) with equation (27) we have 

^ = 47rGao, (31) 

dz 

where 

0-0= Podz. (32) 

Jo 

From equations (25), (28) and (31), we obtain 

1 dPa 



Po dz 



-nj^z - AnGao. (33) 



Linearizeing equations (21) to (26) and using equations (30), (31), (33), we have 

dui dvi dwi 

dx dy dz ' 
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^ + (.0-2^^.^)^ + ^^1 = 0, 



(35) 



dui 



^ + (.0 - 2™--)^ = --^ + ^Cn^-P, + 2C1],., - ^, (36) 

^ + (-0 - |™..)|i + f^... = -^^ - (2C - Im^u, - fi, (37) 

at 2 ay a^; po ay 2 ay 

<9wi ,_ 3^^ ^dwi 1 9Pi 1 9Po <9^i , ^ 

ot 2 ay Po a^; pg oz oz 

dx^ dy"^ dz"^ 
where a perturbed quantity is denoted by subscript 1. 



4. The Effect of the Coriohs Force 

4.1. Eigenvalue Problem 

First, we consider the case without the self-gravity and tidal force, namely, G = 0, 
\t'i = and T = in order to investigate the effect of the Coriolis force. Then we can 
assume that a perturbed quantity /i has the form as 

fi{x,y,z,t) = fi{z) exp[i{k^x + kyy - ujt)], (40) 

where u is the complex angular frequency of perturbed quantities, k^ is the radial wave 
number, and ky is the azimuthal wave number. Note that fi{z) is a function of z. By 
assuming that all the perturbed quantities have the form as equation (40), and the self- 
gravity and the tidal force are neglected, equations (29) and (34)-(38) are rewritten using 
^1 = 0, T = and equations (33) with G = and (40) (we omit "in the following equations) 

vo = -pgV^K/po, (41) 

ikxUi + ikyVi H — ; — = 0, (42) 

dz 

-iuopi + ——wi = 0, (43) 

dz 

-iujui = -ik^—Pi + 2CVLk—Pi + 2CVtKVi, (44) 

Po Po 

—iojvi = —iky — Pi — 'i~u]i ^ 2CQkUi, (45) 

Po dz 
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—iuwi 


1 rfPi fii^ 

, Pi, 
Po dz Po 


uj = uj — kyVQ^z). 


1 dpo 
Pi = — -i-wi. 

loj dz 


e 

dvQ vo dpo 


dz Po dz ' 


d'^vo 
dz^ 


^ / 1 rfpo y 1 rfVo 
Vpo c?2 y Po dz^ 



where 

From Eq. (43), we have 

From equation (41), we have 
and 



From equations (44), (45), (48) and (49), we have 

kxU + 2iCkyQK Pi 



Ml 



u^-AC^Ql Po' 



and 



Vl 



{kyUj — 2iCkx^K 



Pi 

Po 



._dvQ AC^Qj^vodpo 

-iu;-— H — ) wi 

dz luj Po dz 



2n2 \-l 



{f - 4C2l]|) 



Substituting equations (51) and (52) into equation (42), we have 



Pi 



«Po 



{kl + A;2)2c. 



dw\ 



_ dvo AC^Vt]^ vq dpt 



iuj^ - AC^nl)^ + ky { u^ + —^--r- I wi 



(46) 
(47) 
(48) 
(49) 

(50) 



(51) 



(52) 



dz \ dz ijj Po dz 

Substituting equation (48) and (53) into (46), and using equations (49) and (50), we 

d'^wi „dwi 



(53) 



get 



, 2 + E—^ + Fwi = 0, 

dz"^ dz 



where 



E 



1 dpo 
Po dz 



1 + 



sc'n'kyvo 



:io^ 



AC^n^)Lj 



(54) 



(55) 



kyVo 


Vl dpo\ 
VPo dz ) 

1 


Q^- 


- 4C21]2^ 



i_(fpo_ 
Po dz"^ 



{ki+kir(u^+ 



1 dp, 
Po dz 



0o2 



n'z + 



CflxkyVo dpo 
uopo dz 



(56) 
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Thus, perturbation equations (34)-(38) are reduced to a differential equation (54) for wi 

Only odd solutions for Wi are considered because even ones are always stable according 
to our calculations. Thus one of the boundary conditions is given by 

wi = at z = 0. (57) 

Outside the dust layer, from equation (54) to (56), we have 

where 



- K^wi = 0, (5^ 



^x + S/ 



U 



2 



^' - hh^kF,- (5s) 



We select a root K whose real part is positive. Then the outer boundary condition, i.e., 
Wi — s> for z — ^ oo, is satisfied by the solution. 



From equation (60), we have 



Wi oc exp(— i^^;). (60) 



-^ + Kwi = for z > Zd. (61) 

dz 

From equations (53) and (61), we get 

{kl + kl)u dz {kl + kl)u '■ ^'^ 

At the boundary between the dust and the gas layers. Pi and Wi must be continuous. Thus 
equations (53) and (62) read by using equation (49) 

— h A - — — ] wi = at z = Zd. (63) 

dz \ ujpQ dz J 

We obtain eigenvalues u by solving differential equation (54) to satisfy the boundary condi- 
tions (57) and (63). The growth rate of the shear instability is ui = Q{uj), where 0= denotes 
the imaginary part. 

In this work, we employ CRNDDD unlike Ishitsu & Sekiya (2002) which used HDDD. 
However, CRNDDD is not an explicit function of z (see eq. [5]). Thus, we integrate the 
differential equation by using the dust surface density cq as the independent variable. Then 



-12- 



we have to calculate the density po and the height z as functions of the surface density oq. 
From equations (6) and (32) with Po = Pd + Pg, we have 

r 1 dz 

Differentiating equation (5) with respect to u and substituting in equation (64) and inte- 
grating, we get 

(To = v^vrPg In {u/um) + \/{ulu„,f - 1 . (65) 

This equation is solved for u as 

u = Um cosh ( ao/vJrjrpg ) . (66) 

Next, we obtain po using equation (6), 

Po = Pg[um cosh{ao /\/jr]rpg) - q]~^ (67) 

We can get z from equations (5) using equation (66), and Vq from equation (41) using 
equation (67). Differentiating equation (67) by z and using equation (32), we get 

dpo 3 u„ 



Po 



dz \Ur]rp- 



and 



sinh(ao/vO'?7rpg), (68) 



dVo 3 (dpQ\ plum .(in \ (an\ 

-d^ = 7X-dI) - j^^vvf ^^^^"°/"^'''^^- ^''^ 

Thus the differential equation (54) and the boundary condition (63) are rewritten 

d^wi 1 fdpo \ dwi F 

-r^ + — -1- + Po^ H— + —wi = 0, 70 

dao^ Po \dz J dao p^ 

and 

dwi If kyVodpo\ 

-\ \ K - — 1— ] wi = at z = Zd- (71) 



dao Po \ ujpo dz 

If there are an eigenvalue u and an eigenfunction wi of the differential equation (70) under 
the boundary conditions (57) and (71), then their complex conjugates u* and wl also satisfy 
these equations as seen by taking complex conjugates of these equations. Thus, if there is a 
growing mode with the growth rate uj, then there is also a decaying mode with —ui, and 
vice versa. Replacing ky by —ky and u by —u in equations (70) and (71), we get the identical 
equations. Thus, we have a same growth rate of growing modes for opposite signs of ky with 
a same absolute value. Further it is easily seen that eigenvalues are same for opposite signs 
of kx with a same absolute value. Consequently, the growth rate of an unstable mode is an 
even function of k^ and ky. 
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4.2. Results 

The Richardson number J and the dust density at the midplane pd(0) were taken as 
independent parameters in Sekiya & Ishitsu (2000) on determining CRNDDD. However, the 
Richardson number J is determined if the dust density at the midplane Pd{^) and the dust 
surface density factor f^ are given for CRNDDD. Figure 1 draws the Richardson number 
as functions of the dust density at the midplane pd(0) with constant dust surface densities, 
which are calculated using equation (10). 

Figure 2 shows the growth rate uji as a function of the absolute values of radial and 
azimuthal wave numbers kx and ky in the case where Pd{0)/Pg = 1-0, and fa = 1.0 at 
R =1AU. Hereafter, non-dimensional wave numbers, kxfjR and kyfjR, are used instead of the 
dimensional ones, k^ and ky, where ?7i?=constant in our local Cartesian coordinate system. 
This figure shows that there exists a maximum radial wave number for the shear instability 
kxc{> 0) for each value of ky] ii\kx\ > k^c, the growth rate ui = 0. Figure 3 shows the growth 
rate uj as a function of \kx\rjR with keeping \ky\riR =10.5 in the case where PdiS^)l Pg = 1-0, 
and fd = 1.0 (same as Fig. 2). We have kxcf]R =190 from this figure. Figure 4 shows a 
graph similar to Figure 3, but for pd{0)/pg = 40 and \ky\riR =180; in this case, we have 
kxcf]R = 340. These figures will help us to understand the stabilization and amplification 
mechanisms of the shear instability by the tidal force stated in §5. 

Figure 5 shows the growth rate of the mode with the most unstable wave number 
(hereafter referred to as "the peak growth rate") as a function of the dust/gas density ratio 
on the midplane. In the case without the Coriolis force, the peak growth rate is of the order 
of the Kepler angular velocity Qk over the entire dust density and increases monotonically 
for density smaller than Pd{0)/pg = 50 except for 2 < Pd{0)/pg < 10. On the other hand, the 
peak growth rate decreases for pd{0)/Pg > 50. In the case with the Coriolis force, the peak 
growth rate is smaller than that without the Coriolis force. Moreover, for Pd{0)/Pg > 100, 
the shear instability is stabilized! 

This is explained as follows. For a system without the Coriolis force, the energy of 
instability is supplied from the unperturbed flow at the co- rotation sheet, where the angular 
phase velocity of the perturbed quantities is coincident with the unperturbed azimuthal 
velocity (see Sekiya & Ishitsu 2000 for detail). However, when the Coriolis force is included, 
other resonances like the Lindblad resonances appear. In the system with the Coriolis force, 
perturbed energy is supplied from the latter "resonances" by the baroclininic instability 
instead of the co-rotation resonance. Figure 6 shows that the real and imaginary parts of 
the eigenfunction ^{wi) and Q{wi), respectively, in the case where kyrjR = 32.4, kx = 0, 
Pd{0)/pg = 12, C = 1 and fd = 1. The positions of the co-rotation resonance, and the upper 
and the lower resonances are presented in Figure 6. Figure 7 also shows the eigenfunctions in 
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the case with a higher density where kyTjR = 180, k^ = 0, Pd{0)/pg = 40, C = 1 and fd = 1. 
Compared to Figure 6, the intervals between the resonances become narrower as the dust 
density at midplane increases. The amphtude of an eigenfunction largely changes around 
the co-rotation, the upper and the lower resonances. It becomes more difficult to connect the 
eigenfunction as Pd{0)/pg increases because the intervals between the resonances decrease. 
Thus, no eigenvalues exist for Pd{0)/pg > 100, namely the shear instability is stabilized. 



5. The Effects of the Tidal Force and the Self-Gravity 

In this section the linear perturbation equations are solved taking both the tidal force 
and the self-gravity into account as well as the Coriolis force. The tidal force plays a crucial 
role for the stabilization of the shear instability. 



5.1. Formulation 

The normal mode analysis (i.e. the Fourier transform) with respect to x cannot be done 
in the coordinate system written in the previous section because the coefficients of equations 
depend on x, ii T ^ 0. Thus we transform the coordinate y into the coordinate y' shearing 
with the local Kepler velocity —^TQ^x (see Eq. [19]) in order that the normal mode analysis 
for X can be done, 

y' = y + ^TnKxt. (72) 

We assume that perturbed quantities are written 

fi{x,y,z,t) = fi{z,t)exp[i{kyy' + k^x)]. (73) 

Substituting this expression into equations (34)-(39), and abbreviating fi{z,t) as /i, the 
perturbation equations are rewritten 

ik'^ui + ikyVi + —— = 0, (74) 

— + tkyVopi + -1-^1 = 0, (75) 

^' + ikyvoui = -ik'^—Pi + 2CnK—pi + 2CnKVi - ik'^^i, (76) 



dVl ., _ -7 1 n ^^0 { c^^ 3 



ikyVoVi = —iky — Pi — 1""^! ~ ( 2C — -T 1 VLkUi — iky'^i, (77) 

Po dz V 2 / 



dwi , I dPi I dPo 
dt po dz pI dz '^ 


dz 


^ - {k'^ + kD^i = ^-r^Gp,, 




K = K + ^TkynKt. 
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(78) 

(79) 
uz- 

where, 

The coefficients of these equations depend on t. Thus, the normal mode analysis (i.e. the 
Fourier transform with respect to t) cannot be performed in the case with the tidal force. 
Additionally, the self-gravity brings out an extra free parameter when obtaining an eigen- 
value, so that the calculation becomes difficult. Thus, we numerically integrate equations 
(74) to (80) following the method of Ishitsu & Sekiya (2003) (see Appendix A). 

It must be noted that some unperturbed quantities such as po and vq are not differ- 
entiable at Zd for CRNDDD (see Figs. 1 to 3 in Sekiya & Ishitsu 2000). Accordingly our 
equations become singular there. In order to avoid it, we soften density distribution around 
z = Zd a.s follows. We choose a point Zc which satisfies <^ Zc < Zd- For z < Zc, we use 
CRNDDD; on the other hand, for z > Zc, dust density is assumed 

Pd{z) = Aexp[-B{z-z,)], (81) 

where constants A and B are determined so that the dust density is continuous and differ- 
entiable at z = z^. In present paper, we put Zc = 0.95zd in all the calculations. Here we 
adopted this value of Zc so that the value dose not affect the results and our calculations do 
well. Actually we conffimed that our code worked well when Zc = 0.95. Then the Richardson 
number for z > Zc is always larger than the constant Richardson number for z < Zc- 

Boundary conditions for the z-direction are given as follows. The mirror symmetry with 
respect to the midplane is assumed: 

wi = at z = 0, (82) 

f)P 

—1 = at z = 0, (83) 

oz 

and 

_^ = at z = 0. (84) 

oz 

The continuity of the pressure at z = Zd, i.e. the boundary between the dust and the gas 
layers is applied in the case without the tidal force (see §4). However, it is difficult to apply 
this condition to the case with the tidal force. Thus, we solve the perturbation equations 
numerically within region [0, zq], where the solid-wall condition is applied at zq for simplicity: 

Wi = at z = Zq. (85) 
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The value of zq is chosen to be large enough in order for an eigenfunction to decay sufficiently 
at a boundary zq. We used zq = 2zd in our calculations, and confirmed that numerical 
solutions for T = with this condition agree well with the results of the normal mode 
analysis in §4. From equations (75), (82) and (85), we have 

— — h ikyVopi = at z = and Zq. (86) 

ot 

If we give an initial condition with pi = at z = and zq, equation (86) implies that pi = 
for t > 0. Hereafter we take only such initial conditions for simplicity. Accordingly boundary 
conditions for pi are given by 

pi = at z = and zq. (87) 

By assuming that pi = for z > zq, equation (79) reads 



dz 
where 



^'*'-*'^*. = 0, 



k' = ^If + tj. (89) 

The solution of (88), which satisfies ^i ^ for z —^ 0, is given by 

^1 oc exp{-k'z). (90) 

From equation (90), we have the boundary condition for ^i at z = Zq: 

—1 = -A;'^i at z = Zq. (91) 

oz 

From (78), (85), (87) and (91), the boundary condition for Pi at z = zq is given by 

dPi 



dz 



k'po^i at z = Zq. (92) 



5.2. The Effect of the Self-Gravity 

Figure 8 shows the results neglecting both the Coriolis force and the tidal force in order 
to see the effects of the self-gravity of the dust layer. The self-gravity of the dust layer does 
not affect the shear instability as much as the Coriolis force and the tidal force. The self- 
gravity decreases the growth rate of the shear instability somehow. It is natural to expect 
that the self-gravity has little effects on the shear instability by noting that the growth rate 
of the gravitational instability is highest for a much smaller radial wave number k^- Thus, 
the self-gravity does not change the results of the shear instability. 
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5.3. The Effect of the Tidal Force 

Here, we mainly investigate the stability of shear flow for pd{0)/pg < 100 for which the 
Coriolis force cannot stabilize the shear instability (see §4). For Pd{0)/pg = 40, Figure 9 
shows the time evolution of the radial, azimuthal and vertical components of the spatially 
averaged kinetic energy of the perturbed flow; the values of the parameters in Figure 9 are 
same as Figure 4, except that T = 1 in Figure 9, while T = in Figure 4. These components 
vibrate, and the growth rate of the shear instability is effectively zero. Thus, the tidal force 
stabilizes the instability. The calculations for several other parameters have been performed, 
and the tidal force always stabilizes the instability. 

This stabilization is explained by the following mechanism (see also Ishitsu & Sekiya 
2003). If the tidal force does not work (T = 0), we can obtain the growth rate of the shear 
instability with the method written in section 4. The growth rate depends on the radial and 
the azimuthal wave numbers k^ and ky. The radial wave number plays a crucial role. There 
exists a critical radial wave number k^c {> ^)] ^i \kx\ > k^^, then the growth rate is zero (see 
Figs. 2 to 4). If the tidal force is included, that is, the radial shearing works, the radial wave 
number changes with time due to the Keplerian shear: 

k', = k, + ^TkytQK- (93) 

li \kx\ < kxc, the perturbed quantities grow at first. As time passes, the value of |A;^| changes, 
and eventually \k'^\ exceeds the critical wave number kxc, then the growth rate becomes zero, 
that is, the shear instability is stabilized. From equation (93), we evaluate the stabilization 
time tg by 

ts = [sgn{ky)kxc -kx]/ I -TkyQx j ■ (94) 

In the case of Figure 4, kxc^jR = 340 and we have ts = 1.3i7^^ for kx = 0. As seen in Figure 
9, the perturbation energy does not grow for t > tg. 

Additionally, the change of the radial wave number with time has not only the stabilizing 
effect but also amplification mechanism in some occasions. Figures 10 and 11 show time 
evolution of perturbation energy and the growth rate, respectively, in the case where kxTjR = 
—190 and kyrjR = 10.5 (same values of the parameters as Figure 3, except that T = 1 in 
Figures 10 and 11, while T = in Figure 3). The growth rate oscillates, but the averaged 
value for each oscillation increases at first {tflx <12, i.e. k'^ <0). After that time, however, 
it begins to decrease (12< tQx <24, i.e. 0< k'^rjR <190), and finally it stops growing (24< 
tQx, i-e. 190< k'^rjR). These changes of the growth rate of the energy can be interpreted by 
considering the change of the radial wave number in Figure 3. If the value of kx increases from 
the initial value kxf]R = —190 with keeping kyfjR = 10.5, the growth rate first increases, and 
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reaches a peak value at kx = 0; after that, the growth rate decreases, and finally, the mode is 
stabilized for k^riR > k^cVR =190. The increase of the growth rate during 0< tQ^ <12 (i.e. 
k'^ <0) in Figure 9 resembles the 'swing amplification' (Goldreich & Lynden-Bell 1965) in the 
sense that the Keplerian shear causes the growth. However, our equations are different from 
Goldreich & Lynden-Bell (1965) in the points that the fiuid is incompressible and the self- 
gravity plays no important roles in our case. Although the growth rate increases for adequate 
ranges of k'^, the stabilization effect would consequently overcome the amplification. In fact, 
we confirmed by a number of numerical simulations that the amplification has little effect 
on our results. 

In the case without the tidal force, the growth rate is of the order of the Kepler angular 
velocity. The instability with such a small growth rate is stabilized by the tidal force (Ishitsu 
& Sekiya 2003). We have shown in this subsection that CRNDDD is stable against the shear 
instability for pd{0)/Pg < 100. Moreover, CRNDDD is stabilized due to the Coriolis force for 
Pd{^)/Pg > 100 as shown in §4. Thus, CRNDDD is stable for all range of the dust density 
at the midplane. 



6. Discussion and Conclusions 

Recently, Ishitsu & Sekiya (2003) presented that the tidal force has the effect to stabilize 
the shear instability. They used the constant dust distribution with sinusoidal transitional 
zones as an unperturbed state, which corresponds to an initial distribution before the onset 
of the shear instability. They showed that the shear instability develops before the onset 
of the gravitational instability because the growth rate of the shear instability becomes so 
large that the Coriolis and the tidal forces can not stabilize it. Thus, the dust layer would 
become turbulent. The dust density distribution changes due to the turbulent diffusion 
and the Richardson number is expected to approach a constant value. In this paper, we 
used the constant Richardson number dust density distribution as an unperturbed state and 
investigated its linear stability. 

First, we showed that the Coriolis force stabilizes the shear instability if the dust density 
on the midplane is larger than a hundred times the gas density. Next, we found that the 
coupling effect of the tidal and the Coriolis forces stabilizes the shear instability also for 
the midplane dust density less than a hundred times the gas density. Thus, the constant 
Richardson number dust density distribution is always stable against the shear instability. 
Dust settling proceeds, and eventually the dust density on the midplane exceeds the critical 
density of the gravitational instability. It is possible that the gravitational instability of 
the dust layer occurs even for the dust to gas surface density ratio is that for the solar 
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composition. So far, it has been considered that the gravitational instabihty requires the 
special mechanisms of the dust concentration or the gas dissipation in order to overcome the 
shear instability (Sekiya 1998; Youdin & Shu 2002; Youdin & Chiang 2004). The stabilizing 
effects of the tidal and the Coriolis forces elucidated in this paper relaxes this restriction. 
And the gravitational instability may be a common process of the planetesimal formation. 

In this paper, we assumed that the constant dust density distribution as an unperturbed 
state, and showed that it is stable. In the real evolution of the dust layer, the dust settling 
and the turbulent mixing may cause the dust distribution somewhat different from that with 
the constant Richardson number. It is also unclear yet whether the planetesimals really form 
after the dust density exceeds the critical density of the gravitational instability. In order 
to elucidate the actual evolution of the dust layer, we plan to perform nonlinear, dust-gas 
multifluid numerical simulation in future. 

We thank Dr. Seiichiro Watanabe for valuable comments. The calculations in this work 
were partly performed with computers at Astronomical Data Analysis Center, National 
Astronomical Observatory of Japan. 



A. Numerical Method 

We adopt the MAC method, in which the pressure is determined by the condition that 
the equation of continuity is satisfied in the next step, and other variables Ui,Vi,Wi and pi 
are calculated using this value of the pressure. 

We define the divergence of the perturbed velocity by 

D = ik'^ui + ikyVi + — — . (Al) 

Multiplying equations (76) by ik'^ and equation (77) by iky, and taking partial derivative of 
equation (78) with respect to z, and adding, we have 

^^ = -tkyVoD-2tky^wi- — [-k'^-kl + —^]Pi + -^-^^ 



dt ' " " dz^ po\ "" y dzy ' pldz dz 

+2iCnKk'^—pi + 2iCnKk'^vi -i{2C- 3T) ^KkyUi 
Po 

+^(\^Pi)-^nGp,. (A2) 



(A3) 



dz \pI dz 
We use the first order approximation: 

dD D"+i - D 



dt At 
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We require that the equation of continuity (74) is satisfied in the next step, i.e. D"'^^ = 0. 
Thus, we get the equation for Pi: 

,n2 ,2 , -92 1 dpQ d 



— i-'"^ — i-^ -I- ^•~r^_~_ pn 



dz^ Po dz dz ^ 
po I ( -^ - ikyvA D'' - 2iky-^w'l + 2iVtKCk'^vl -i{2C- 3T) VtKkyUl 



Jn^O .„ , d f 1 dPo 



+2.Cn,-C^ri- + -^-,^prJ-4.Gp;-|^ (A4) 

From this equation, we can get Pi at the time step n. We calculate ^ at the time step n 
from 

f)2^n 

-^ - {k'f + k^^^i = 47^Gp^ (A5) 

Next, we calculate pi,ui,vi and wi at the time step n + 1 from equations (75) to (78) 
using approximation like (A3): 



p^^ = p^ + At I -^A;,t;op^ " ^"^^ / ' ^^^^ 

M^+i = m;^ + At (-iA;^t;o< - ik'^—Pl' + 2CnK—Pi + 2Cl]^t;^ - i/t^^^M , (A7) 

I Po Po J 



,«+!_„,« I A+J .•/.„-,„,« 0_^,„ ., 1 p„ I ^^ 3, 



t;5^+^ = t;^ + At <^ -iA;j,t;ot;r - -f^< - iky— PI" - ( 2C - -T 1 fixw" - «A;j,^^ ^ , (A8) 

u,« ^ < . A* {-.*,..< - If: + -i^rf. - f } . (A9) 

The above method is the first order accuracy with respect to At. In order to improve them 
to the second order accuracy, we use the following strategy. We replace perturbed quantities 
/" = (p", -u", Vi, Wi) on right hand side of equation (A4) with 

J.n^ ^ ^^„+i ^ ^„^^2. (AlO) 

Again, we solve equation (A4) replacing n by n + 1/2. More exact values of p"^^, Ui'^'^, f "^^ 
and Wi~^^ are obtained by replacing these quantities at n in braces of the right hand sides of 
equations (A6) to (A9) with ones at n + 1/2. We adopt one dimension staggered mesh, where 
Wi is estimated at grids and Pi, ^i,pi,Mi,fi at midpoints of adjacent meshes. In equations 
(A4), (A6) and (AS), wi is calculated by taking the mean values at the adjacent meshes. So 
as pi in equation (A9). 

Initial conditions are set as follows. Some Fourier components of lower orders are selected 
as to satisfy the boundary conditions of Mi, Wi and pi. Each Fourier coefficient is given by 
a random number. Velocity Vi is determined from Ui and Wi by using the equation of 
continuity (74). 
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Fig. 1. — The Richardson number J as functions of the dust /gas density ratio on the mid- 
plane pd{0)/pg for fd = i (sold hue) and f^ = 17.15 (dotted line) at R =1AU. The position 
of the arrow shows the critical dust density of the gravitational instability. The dot-dashed 
line denotes the critical Richardson number X. = 0.25. 
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Fig. 2. — The growth rate uj of the mode as a function of the absolute values of the radial 
and the azimuthal wave numbers \kx\ and \ky\, respectively, in the case where C = 1, T = 
and G = with J = 4.55 x 10"^ Pd{0)/pg = 1.0 and fa = 1.0 at R =1AU. 
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Fig. 3. — The growth rate ujj as functions of \kx\ in the case where \ky\riR = 10.5 with 
J = 4.55 X 10-2, pdiO)/pg = I, C = I, T = and G = 0. 
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Fig. 4. — The growth rate ujj as functions of \kx\ in the case where \ky\r]R = 180 with 
J = 1.18 X 10-3, Pd{0)/pg = 40, C = 1, T = and G = 0. 
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Fig. 5. — The growth rate ujj of the mode with the most unstable wave number as functions 
of the dust/gas density ratio on the midplane pd{0)/pg, with the Coriohs force (sohd curve) 
and without this force (dotted curve). 
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Fig. 6. — The real and the imaginary parts of the eigenfunction ^{wi) and '<s{wi), respec- 
tively, in the case where kyTjR = 32.4, k^ = 0, pd{0)/pg = 12, C = 1, T = 0, G = and 
/d = 1 at i? =1AU. The horizontal line denoted by CR shows the co- rotation sheet. The 
horizontal lines denoted by UR and LR show the upper and lower resonances, respectively. 



-29- 



3[Wi]/77V^ 



0.03- 



^ 
> 



0.02- 



0.01- 







k^r]R=\m 



K=o 



/=L18xlO 
R=IAV 

c=\ 

T=0 
G=0 



-3 



-5 




5 



UR 
CR 
LR 



Fig. 7. — The real and imaginary parts of the eigenfunction ^{wi) and '<s{wi), respectively, 



in the case where kyrjR 
R =1AU. 



180, k^ = 0, pd{0)/pg = 40, C = 1, T = 0, G = and /rf = 1 at 
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Fig. 8. — The time evolution of the growth rate in the case including {G ^ 0, sold line) 
and neglecting (G = 0, dotted line) the self-gravity with J = 9.48 x 10"'', Pd{0)/pg = 100, 



C = 0, T = 0, and kyrjR 
in section 4. 



394. The dotted line shows the solution derived by the method 
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Fig. 9. — The time evolution of radial, azimuthal and vertical components of the the spatially 
averaged kinetic energy of the perturbed flow, in the case including the tidal force (T = 1) 
with J = 1.17 X 10-3, pdiO)/pg = AO, C = 1, G ^ 0, K = and kyfiR = 180 at R =1AU. 
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Fig. 10. — The time evolution of radial, azimuthal and vertical components of the the 
spatially averaged kinetic energy of the perturbed flow, in the case including the tidal force 
(T = 1) with J = 4.55 x 10"^ PdiO)/Pg = I, C = I, G ^ 0, Kt]R = -190 and kyfiR = 10.5 
at R =1AU. 
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Fig. 11. — Time evolution of tlie growtli rate in tlie case including tidal force (T — 1) with 
J = 4.55 X 10-^ PdW/pg = 1, C = 1, G ^ 0, ti(fl = -190 and kyi^R = 10.5 at R =1AU. 



